# -*- coding: utf-8 -*-
"""
    @project: pythonProject
    @Author：HanYonghua
    @file： sinc_interpolation3.py
    @date：2025/7/24 19:54
    @blogs: https://www.ncatest.com.cn
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import time

def sinc_resample_enhanced(x, fs_old, fs_new, L=8, window='hamming', anti_alias=True):
    """
    增强版Sinc插值函数
    参数:
        x: 输入信号
        fs_old: 原采样率
        fs_new: 目标采样率
        L: 插值半长
        window: 窗类型 ('rect', 'hamming', 'hann', 'blackman')
        anti_alias: 是否启用抗混叠滤波
    返回:
        重采样后的信号
    """
    # 抗混叠滤波（升采样时必需）
    if anti_alias and fs_new > fs_old:
        cutoff = 0.9 * (fs_old / 2)
        b = signal.firwin(101, cutoff, fs=fs_old)
        x = signal.lfilter(b, 1, x)

    ratio = fs_new / fs_old
    num_samples_new = int(len(x) * ratio)
    t_new = np.arange(num_samples_new) / fs_new
    x_new = np.zeros(num_samples_new)

    # 窗函数选择
    if window == 'hamming':
        win = lambda u: 0.54 + 0.46 * np.cos(np.pi * u / L)
    elif window == 'hann':
        win = lambda u: 0.5 + 0.5 * np.cos(np.pi * u / L)
    elif window == 'blackman':
        win = lambda u: 0.42 + 0.5 * np.cos(np.pi * u / L) + 0.08 * np.cos(2 * np.pi * u / L)
    else:  # rect
        win = lambda u: 1.0

    # 边界安全处理
    pad_width = L
    x_padded = np.pad(x, (pad_width, pad_width), mode='edge')

    for i, t in enumerate(t_new):
        n_center = int(t * fs_old) + pad_width
        delta = (t - (n_center - pad_width) / fs_old) * fs_old
        u = np.arange(-L, L + 1) - delta
        coeffs = np.sinc(u) * win(u)
        x_new[i] = np.sum(x_padded[n_center - L: n_center + L + 1] * coeffs)

    return t_new, x_new


# 测试信号生成
fs_old = 44100
t_old = np.arange(0, 0.02, 1 / fs_old)  # 20ms信号
f_signal = 1000  # 1kHz正弦波 + 5kHz高频
x_old = np.sin(2 * np.pi * f_signal * t_old) + 0.3 * np.sin(2 * np.pi * 5000 * t_old)

# 重采样执行
fs_new = 96000
start_time = time.time()
t_new, x_new = sinc_resample_enhanced(x_old, fs_old, fs_new, L=8, window='blackman')
print(f"计算耗时: {time.time() - start_time:.4f}秒")


# 频谱分析函数
def plot_spectrum(x, fs, label, color):
    f, Pxx = signal.welch(x, fs, nperseg=1024)
    plt.semilogy(f, Pxx, label=label, color=color, linewidth=1.5)


# 设置中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 频谱对比
plt.figure(figsize=(12, 5))
plot_spectrum(x_old, fs_old, '原始信号 (44.1kHz)', 'blue')
plot_spectrum(x_new, fs_new, '插值信号 (96kHz)', 'red')
plt.xlim(0, 20000);
plt.ylim(1e-8, 1)
plt.title("频谱对比（抗混叠滤波+Blackman窗）", fontsize=14)
plt.xlabel("频率 (Hz)");
plt.ylabel("功率谱密度")
plt.legend();
plt.grid(True)

# 时域对比
plt.figure(figsize=(12, 5))
plt.plot(t_old, x_old, 'bo', label='原始信号', markersize=4, alpha=0.5)
plt.plot(t_new, x_new, 'r.', label='插值点', markersize=6, alpha=0.8)
plt.xlim(0.01, 0.012)  # 聚焦10-12ms区间
plt.title("时域信号对比（插值点显示）", fontsize=14)
plt.xlabel("时间 (s)");
plt.ylabel("幅值")
plt.legend();
plt.grid(True)

plt.tight_layout()
plt.show()